/*
  File:    Microenvironment.cpp
  Purpose: Utilities to update the microenvironment, including:
               matrix degrading enzyme
               fibronectin
               structure
               mobility?

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>

using namespace std;

#include "../ImprovedMath.h"
#include "../Matrix.h"
#include "../DiffusionSolver.h"

#include "Polynomial.h"

Polynomial MDEdiffusivityPolynomial;

Matrix Original_ECM; 

bool CreateMDEdiffusivityPolynomial( void )
{
 extern double MDEdiffusionConstant;
 extern double MaxMDEdiffusionConstant;
 
 Matrix Coeff = zeros(4,4);
 Matrix RHS = zeros(1,4);
 
 Matrix Soln = zeros(1,4);
 
 *Coeff(0,0) = 1;
 *Coeff(0,1) = 0.0;
 *Coeff(0,2) = 0.0;
 *Coeff(0,3) = 0.0;
 *RHS(0) = MaxMDEdiffusionConstant;
 
 *Coeff(1,0) = 0;
 *Coeff(1,1) = 1.0;
 *Coeff(1,2) = 0;
 *Coeff(1,3) = 0;
 *RHS(1) = 0;
 
 *Coeff(2,0) = 1.0;
 *Coeff(2,1) = 1.0;
 *Coeff(2,2) = 1.0;
 *Coeff(2,3) = 1.0;
 *RHS(2) = MDEdiffusionConstant;
 
 *Coeff(3,0) = 0.0;
 *Coeff(3,1) = 1.0;
 *Coeff(3,2) = 2.0;
 *Coeff(3,3) = 3.0;
 *RHS(3) = 0.0;
 
 Soln = GaussianEliminationPartialPivoting( Coeff, RHS );
 
 MDEdiffusivityPolynomial.SetDegree(3);
 for( int i=0; i <= MDEdiffusivityPolynomial.GetDegree() ; i++ )
 { MDEdiffusivityPolynomial.SetCoefficient( i, *Soln(i) ); }
 
 return true;
}

bool UpdateMDEdiffusivity( void )
{
 extern int Xnodes;
 extern int Ynodes;

 extern Matrix MDEdiffusivity;
 extern Matrix ECM;
 
 extern double MDEdiffusionConstant;
 
 extern Matrix Structure;

 if( MDEdiffusivity.TellRows() != Xnodes ||
     MDEdiffusivity.TellCols() != Ynodes )
 { MDEdiffusivity = ones(Xnodes,Ynodes); }
 
 for( int i=0 ; i < MDEdiffusivity.TellRows() ; i++ )
 {
  for( int j=0 ; j < MDEdiffusivity.TellCols() ; j++ )
  {
   double ecm = *ECM(i,j);
   if( ecm >= 1 ) 
   { ecm = 1.0; }
   if( ecm < 0.0 )
   { ecm = 0.0; }
   *MDEdiffusivity(i,j) = MDEdiffusivityPolynomial(ecm);
   
   if( *Structure(i,j) > 0.1 )
   { *MDEdiffusivity(i,j) = 0.001; }
   
  }
 }
 return true;
}

bool IterateBVD( double ds )
{
 extern bool Disable_bulk_o2_source_degradation;
 if( Disable_bulk_o2_source_degradation )
 { return true; }
 
 extern double _EPS;

 extern Matrix BloodVesselDensity;
 extern Matrix MDE;
 
 extern double ECM_DegradationConstant;
 extern double Bulk_o2_source_degradation_rate;
 
 for( int i=0 ; i < BloodVesselDensity.TellRows() ; i++ )
 {
  for( int j=0 ; j < BloodVesselDensity.TellCols() ; j++ )
  {
   double bvd = *BloodVesselDensity(i,j);
   double m = *MDE(i,j);
   
//   double bvd_change = -ds*(ECM_DegradationConstant)*(m*bvd);
   double bvd_change = -ds*(Bulk_o2_source_degradation_rate)*(m*bvd);

   *BloodVesselDensity(i,j) = bvd + bvd_change; // - ds*(ECM_DegradationConstant)*(m*ecm);
   if( *BloodVesselDensity(i,j) < 0.0 )
   { *BloodVesselDensity(i,j) = 0.0; }
   
  }
 }
 return true;
}

bool IterateECM( double ds )
{
 extern double _EPS;

 extern Matrix ECM;
 extern Matrix MDE;
 extern Matrix MDEdiffusivity;
 
 extern double ECM_DegradationConstant;
 
 extern double ECM_ProductionConstant;
 
 extern Matrix Phi, PhiNecrotic;
 extern double _EPS;
 

 
 for( int i=0 ; i < MDEdiffusivity.TellRows() ; i++ )
 {
  for( int j=0 ; j < MDEdiffusivity.TellCols() ; j++ )
  {
   double ecm = *ECM(i,j);
   double m = *MDE(i,j);
   
   double ecm_change = -ds*(ECM_DegradationConstant)*(m*ecm);

   if( *Phi(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS )
   { ecm_change += ds*ECM_ProductionConstant*(1-ecm); }
   
   *ECM(i,j) = ecm + ecm_change; // - ds*(ECM_DegradationConstant)*(m*ecm);
   if( *ECM(i,j) < 0.0 )
   { *ECM(i,j) = 0.0; }
   
   double original_ecm = *Original_ECM(i,j);
   double original_ecm_change = -ds*(ECM_DegradationConstant)*(m*original_ecm);
   
   *Original_ECM(i,j) = original_ecm + original_ecm_change; 
   if( *Original_ECM(i,j) < 0.0 )
   //if(*Phi(i,j)<= _EPS)
   { *Original_ECM(i,j) = 0.0; }

  }
 }
 return true;
}



double MDEsourceFunction(double mde, double x, double y,int i,int j)
{
 extern double _EPS;

 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix ECM;
 
 extern double MDEprodctionConstant;
 extern double MDEdecayConstant;
 extern double MDE_DegradationConstant;

 double output = 0.0;
 
 // no interpolation is needed because source is only called at node points
 
 double ecm = *ECM(i,j);

 if( *Phi(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS )
 { output += MDEprodctionConstant*(1.0-mde); } // MDEprodctionConstant; }
 output -= MDEdecayConstant*mde;
 output -= mde*ecm*MDE_DegradationConstant;

 return output;
}

double MDEdiffusionCoefficientFunction( double mde,double x, double y,int i,int j)
{
 extern Matrix MDEdiffusivity;
 
 extern Matrix X;
 extern Matrix Y;

 double theta_x = x - *X(i);
 if( theta_x < 0 )
 { i--; theta_x = x - *X(i); }
 double theta_y = y - *Y(j);
 if( theta_y < 0 )
 { j--; theta_y = y - *Y(j); }
 
 return (1.0-theta_x)*( (1.0-theta_y)*(*MDEdiffusivity(i,j    )) 
                       +(    theta_y)*(*MDEdiffusivity(i,j+1  )) )
       +(    theta_x)*( (1.0-theta_y)*(*MDEdiffusivity(i+1,j  )) 
                       +(    theta_y)*(*MDEdiffusivity(i+1,j+1)) );
}

double ZeroFunctionBC_MDE( double )
{ return 0.0; }





bool IterateMDE( double ds , double Tolerance )
{
 extern double _EPS;

 extern Matrix X;
 extern Matrix Y;
 extern int Xnodes;
 extern int Ynodes;
 
 extern Matrix MDE;

 MDE = NonlinearDiffusionIteration( 
  X,Y,
  MDE,ds, 
  ZeroFunctionBC_MDE, 'N',  
  ZeroFunctionBC_MDE, 'N',     
  ZeroFunctionBC_MDE, 'N',   
  ZeroFunctionBC_MDE, 'N',  
  MDEdiffusionCoefficientFunction, 
  MDEsourceFunction,
  10000, Tolerance);
  
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0 ; j < Ynodes ; j++ )
  {
/*  
   if( *MDE(i,j) <= -_EPS )
   {
    *MDE(i,j) = 0.0;  
   }
*/   
   if( *MDE(i,j) < 0.0 )
   { *MDE(i,j) = 0.0; }  
  }
 } 
  
 return true;  
}



double DRGsourceFunction(double DRG, double x, double y,int i,int j)
{
 extern double _EPS;

 extern Matrix BloodVesselIndicator;
 extern Matrix TracerT;
 extern Matrix Phi, PhiNecrotic;
 extern Matrix Oxygen;
 
 extern double DRGprodctionConstant;
 extern double DRGdecayConstant;
 extern double DRG_DegradationConstant;

 double output = 0.0;
 
 // no interpolation is needed because source is only called at node points
 //if(*TracerT(i,j)>DRG)
 if(*TracerT(i,j)>0)
 {
 output += DRGprodctionConstant;
 //output += DRGprodctionConstant*(*TracerT(i,j)-DRG);
 // 
 
 }
 
 
 output -= DRGdecayConstant*DRG;
 
 
 if( *Phi(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS )
 { 
 output -= DRG*DRG_DegradationConstant;
 } // MDEprodctionConstant; }
 
 

 return output;
}

double DRGsourceFunctionP(double DRG, double x, double y,int i,int j)
{
 extern double _EPS;

 extern Matrix BloodVesselIndicator;
 extern Matrix TracerS;
 extern Matrix Phi, PhiNecrotic;
 extern Matrix Oxygen;
 extern double TracerMax;
 extern Matrix PreExist;
 extern double DRGprodctionConstant;
 extern double DRGdecayConstant;
 extern double DRG_DegradationConstant;

 double output = 0.0;
 
 // no interpolation is needed because source is only called at node points
 //if(*TracerT(i,j)>DRG)
 if(*TracerS(i,j)>0&&(*PreExist(i,j)==0)&&(*BloodVesselIndicator(i,j)>0)&&(TracerMax>0))
 {
 //output += DRGprodctionConstant;
 output += DRGprodctionConstant*(*TracerS(i,j)-DRG);
 // 
 
 }
 
 
 output -= DRGdecayConstant*DRG;
 
 
 if( *Phi(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS )
 { 
 output -= DRG*DRG_DegradationConstant;
 } // MDEprodctionConstant; }
 
 

 return output;
}


double DRGdiffusionCoefficientFunction( double DRG,double x, double y,int i,int j)
{
extern double D_DRG;
return D_DRG;
}

double ZeroFunctionBC_DRG( double )
{ return 0.0; }

bool IterateDRG_M( double ds , double Tolerance )
{
 extern double _EPS;

 extern Matrix X;
 extern Matrix Y;
 extern int Xnodes;
 extern int Ynodes;
 extern double DRGprodctionConstant;
 extern Matrix TracerD;
 extern Matrix TracerT;
 // cout<<"start"<<endl;
 // getchar();

			
 TracerD = NonlinearDiffusionIteration( 
  X,Y,
  TracerD,ds, 
  ZeroFunctionBC_DRG, 'N',  
  ZeroFunctionBC_DRG, 'N',     
  ZeroFunctionBC_DRG, 'N',   
  ZeroFunctionBC_DRG, 'N',  
  DRGdiffusionCoefficientFunction, 
  DRGsourceFunction,
  10000, Tolerance);
  	for( int i=0 ; i < Xnodes ; i++ )
	{
		for( int j=0 ; j < Ynodes ; j++ )
		{
			//if(*TracerT(i,j)>*TracerD(i,j))
			if(*TracerT(i,j)>0)
			{
				
				
				//*TracerT(i,j) -= DRGprodctionConstant*(*TracerT(i,j)-*TracerD(i,j))*ds; 
				
				*TracerT(i,j) -= DRGprodctionConstant*ds;
			}
			else
			{
				*TracerT(i,j) = 0.0;
			}
			
			if( *TracerT(i,j) < 0.0 )
			{ *TracerT(i,j) = 0.0; } 
			
			if( *TracerD(i,j) < 0.0 )
			{ *TracerD(i,j) = 0.0; }  
			
		}
	} 
	  
 // cout<<ds<<endl;
 // getchar();
  



  
 return true;  
}

bool IterateDRG_P( double ds , double Tolerance )
{
 extern double _EPS;

 extern Matrix X;
 extern Matrix Y;
 extern int Xnodes;
 extern int Ynodes;
 double TOL=0.1;
 extern double TracerMax;
 extern double DRG_Effect;
 extern double DRGDecayVessel;
 extern double DRGprodctionConstant;
 extern Matrix TracerD;
 extern Matrix TracerS;
 if(Tolerance>(TOL/(DRG_Effect*TracerMax)))
 {
 Tolerance = (TOL/(DRG_Effect*TracerMax))*0.01;
 }
 // cout<<"start"<<endl;
 // getchar();

for( int i=0 ; i < Xnodes ; i++ )
	{
		for( int j=0 ; j < Ynodes ; j++ )
		{
			if(*TracerD(i,j)<(TOL/(DRG_Effect*TracerMax)))//YOU MAY SCALE THIS TO TOLERANCE!
			//if(*TracerS(i,j)>0)
			
			{
				
				*TracerD(i,j) = 0;
				
				
			}
						
			 
			
		}
	} 
			
 TracerD = NonlinearDiffusionIteration( 
  X,Y,
  TracerD,ds, 
  ZeroFunctionBC_DRG, 'N',  
  ZeroFunctionBC_DRG, 'N',     
  ZeroFunctionBC_DRG, 'N',   
  ZeroFunctionBC_DRG, 'N',  
  DRGdiffusionCoefficientFunction, 
  DRGsourceFunctionP,
  10000, Tolerance);
  
  	for( int i=0 ; i < Xnodes ; i++ )
	{
		for( int j=0 ; j < Ynodes ; j++ )
		{
			if(*TracerS(i,j)>*TracerD(i,j))//YOU MAY SCALE THIS TO TOLERANCE!
			//if(*TracerS(i,j)>0)
			
			{
				
				
				
				
				*TracerS(i,j) -= (DRGprodctionConstant*(*TracerS(i,j)-*TracerD(i,j)) + DRGDecayVessel*(*TracerS(i,j)))*ds; 
				//*TracerS(i,j) -= DRGprodctionConstant*(*TracerS(i,j)-*TracerD(i,j))*ds; 
			}
			else
			{
				*TracerS(i,j) = 0.0;
			}
			
			if( *TracerS(i,j) < 0.0 )
			{ *TracerS(i,j) = 0.0; } 
			
			if( *TracerD(i,j) < 0.0 )
			{ *TracerD(i,j) = 0.0; }  
			
		}
	} 
	
	
	
	  
 // cout<<ds<<endl;
 // getchar();
  



  
 return true;  
}


bool UpdateMicroenvironment( void )
{
 extern double dt;
 extern double dx;
 
 extern double MDEdiffusionConstant;
 extern double MDEprodctionConstant;
 
 extern double MDEdecayConstant;
 extern double MDE_DegradationConstant;
 

 double s = 0;

 extern double ECM_DegradationConstant;

 if( ECM_DegradationConstant < 1e-8 )
 {
  cout << "Fibronectin degradation disabled; MDE and F updates are pointless." << endl;
  return false;
 }
 
 // determine a good ds
 
 int factor = 1;
 double ds = dt/factor;
 
 double extra_consideration = MAX( MDEprodctionConstant, MDEdecayConstant );
 extra_consideration = MAX( extra_consideration , MDE_DegradationConstant );
 
 extra_consideration /= 10.0; 
 if( extra_consideration < 1.0 )
 { extra_consideration = 1.0; }
 
 while( ds > square(dx) )
 { factor++; ds = dt/factor; }

 // determine the stopping criterion
 
 double Tolerance = IntPow(dx,6); 
/* 
 double scale = MIN(MDEdiffusionConstant,MDEprodctionConstant);
 
 if( scale < 1 )
 { Tolerance *= scale; }
*/ 
 
 for( int count = 0 ; count < factor ; count++ )
 {
  cout << count << " of " << factor-1 << endl;
  // update MDE diffusivity 
  
  UpdateMDEdiffusivity();
 
  // update MDE
  
  IterateMDE( ds , Tolerance );
  
  //IterateDRG_P( ds , Tolerance );
  
  // update fibronectin 
  
  IterateECM( ds );
  
  // degrade pre-existing bulk vasculature
  
  IterateBVD( ds );
  
  s += ds;
 } 
 return true;
}
